DESeq2 analyses steps

This section uses the counts data generated in Section 1 to differential expression analyses using DESeq2.

Prerequisites

Packages required for this analyses are loaded as follows:

setwd("~/github/amnion.vs.other_RNASeq")
# load the modules
library(sva)
library(tidyverse)
library(DESeq2)
library(vsn)
library(pheatmap)
library(ggrepel)
library(RColorBrewer)
library(reshape2)
require(biomaRt)
library(EnhancedVolcano)
library(TissueEnrich)
library(plotly)
library(DT)
library(cowplot)
library(biomaRt)

Import datasets

The counts data and its associated metadata (coldata) are imported for analyses

counts = '~/github/amnion.vs.other_RNASeq/assets/counts-subset-v4.txt'
groupFile = '~/github/amnion.vs.other_RNASeq/assets/batch-subset-v4.txt'
coldata <- read.csv(groupFile, row.names=1, sep="\t", stringsAsFactors = TRUE)
cts <- as.matrix(read.csv(counts,sep="\t",row.names="gene.ids"))

visualize the coldata

DT::datatable(coldata)

order and inspect if all the coldata matches the headers of the counts data.

colnames(cts)
#>  [1] "Naive_H9_hESCs_1"                "Naive_H9_hESCs_2"               
#>  [3] "nTE_D1.Naive_H9_hESCs_1"         "nTE_D1.Naive_H9_hESCs_2"        
#>  [5] "nTE_D2.Naive_H9_hESCs_1"         "nTE_D2.Naive_H9_hESCs_2"        
#>  [7] "nTE_D3.Naive_H9_hESCs_1"         "nTE_D3.Naive_H9_hESCs_2"        
#>  [9] "nCT_P3.Naive_H9_hESCs_1"         "nCT_P3.Naive_H9_hESCs_2"        
#> [11] "nCT_P10.Naive_H9_hESCs_1"        "nCT_P10.Naive_H9_hESCs_2"       
#> [13] "nCT_P15.Naive_H9_hESCs_1"        "nCT_P15.Naive_H9_hESCs_2"       
#> [15] "cR_nCT_P15.Naive_H9_hESCs_1"     "cR_nCT_P15.Naive_H9_hESCs_2"    
#> [17] "nCT_P15.409B2_iPSC_hESCs_1"      "nCT_P15.409B2_iPSC_hESCs_2"     
#> [19] "Placenta.derived_tbSCs_CT30_Ex1" "Placenta.derived_tbSCs_CT30_Ex2"
#> [21] "nST.Naive_H9_Ex1"                "nST.Naive_H9_Ex2"               
#> [23] "nEVT.Naive_H9_Ex1"               "nEVT.Naive_H9_Ex2"              
#> [25] "Primed_H9_hESCs_1"               "Primed_H9_hESCs_2"              
#> [27] "pBAP_D1.Primed_H9_hESCs_1"       "pBAP_D1.Primed_H9_hESCs_2"      
#> [29] "pBAP_D2.Primed_H9_hESCs_1"       "pBAP_D2.Primed_H9_hESCs_2"      
#> [31] "pBAP_D3.Primed_H9_hESCs_1"       "pBAP_D3.Primed_H9_hESCs_2"      
#> [33] "CytoTB_7_gestational_wks_1"      "CytoTB_7_gestational_wks_2"     
#> [35] "CytoTB_9_gestational_wks_1"      "CytoTB_11_gestational_wks_1"    
#> [37] "hESC_H1_STB_gt70um_D8_BAP_1"     "hESC_H1_STB_gt70um_D8_BAP_2"    
#> [39] "hESC_H1_STB_gt70um_D8_BAP_3"     "hESC_H1_STB_40.70um_D8_BAP_1"   
#> [41] "hESC_H1_STB_40.70um_D8_BAP_2"    "hESC_H1_STB_40.70um_D8_BAP_3"   
#> [43] "hESC_H1_STB_lt40um_D8_BAP_1"     "hESC_H1_STB_lt40um_D8_BAP_2"    
#> [45] "hESC_H1_STB_lt40um_D8_BAP_3"     "hESC_H1_D8_MEF.CM.and.FGF2_1"   
#> [47] "hESC_H1_D8_MEF.CM.and.FGF2_2"    "hESC_H1_D8_MEF.CM.and.FGF2_3"
all(rownames(coldata) %in% colnames(cts))
#> [1] TRUE
cts <- cts[, rownames(coldata)]

Batch correction

We will use the combat seq (SVA package) to run batch correction based on the bioproject (dataset origin)

cov1 <- as.factor(coldata$BioProject)
adjusted_counts <- ComBat_seq(cts, batch=cov1, group=NULL)
#> Found 2 batches
#> Using null model in ComBat-seq.
#> Adjusting for 0 covariate(s) or covariate level(s)
#> Estimating dispersions
#> Fitting the GLM model
#> Shrinkage off - using GLM estimates for parameters
#> Adjusting the data
all(rownames(coldata) %in% colnames(cts))
#> [1] TRUE
cts <- cts[, rownames(coldata)]

run DESeq2

The batch corrected read counts are then used for running DESeq2 analyses

dds <- DESeqDataSetFromMatrix(countData = adjusted_counts,
                              colData = coldata,
                              design = ~ condition)
vsd <- vst(dds, blind=FALSE)
keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]
dds <- DESeq(dds)
dds
#> class: DESeqDataSet 
#> dim: 16496 48 
#> metadata(1): version
#> assays(4): counts mu H cooks
#> rownames(16496): ENSG00000000003.15 ENSG00000000005.6 ...
#>   ENSG00000288695.1 ENSG00000288698.1
#> rowData names(106): baseMean baseVar ... deviance maxCooks
#> colnames(48): cR_nCT_P15.Naive_H9_hESCs_1 cR_nCT_P15.Naive_H9_hESCs_2
#>   ... hESC_H1_STB_lt40um_D8_BAP_2 hESC_H1_STB_lt40um_D8_BAP_3
#> colData names(3): BioProject condition sizeFactor

Various contrasts are setup as follows (a total of 9 combinations)

res.nCTvsBAP <- results(dds, contrast=c("condition",
                                        "nCT_P10.Naive_H9_hESCs",
                                        "pBAP_D3.Primed_H9_hESCs"))
res.nTEvsSTB <- results(dds, contrast=c("condition",
                                        "nTE_D3.Naive_H9_hESCs",
                                        "hESC_H1_STB_gt70um_D8_BAP"))
res.nCTvsSTB <- results(dds, contrast=c("condition",
                                        "nCT_P10.Naive_H9_hESCs",
                                        "hESC_H1_STB_gt70um_D8_BAP"))
res.nTEvsBAP <- results(dds, contrast=c("condition",
                                        "nTE_D3.Naive_H9_hESCs",
                                        "pBAP_D3.Primed_H9_hESCs"))
res.nTEvsL40 <- results(dds, contrast=c("condition",
                                        "nTE_D3.Naive_H9_hESCs",
                                        "hESC_H1_STB_lt40um_D8_BAP"))
res.nCTvsL40 <- results(dds, contrast=c("condition",
                                        "nCT_P10.Naive_H9_hESCs",
                                        "hESC_H1_STB_lt40um_D8_BAP"))
res.BAPvsL40 <- results(dds, contrast=c("condition",
                                        "pBAP_D3.Primed_H9_hESCs",
                                        "hESC_H1_STB_lt40um_D8_BAP"))
res.BAPvsSTB <- results(dds, contrast=c("condition",
                                        "pBAP_D3.Primed_H9_hESCs",
                                        "hESC_H1_STB_gt70um_D8_BAP"))
res.STBvsL40 <- results(dds, contrast=c("condition",
                                        "hESC_H1_STB_gt70um_D8_BAP",
                                        "hESC_H1_STB_lt40um_D8_BAP"))

function to save results

processDE <- function(restable, fileName) {
  restable <- restable[order(restable$padj), ]
  outTable <- merge(as.data.frame(restable), 
                    as.data.frame(counts(dds, normalized=TRUE)), 
                    by="row.names", sort=FALSE)
  names(outTable)[1] <- "Gene"
  newName=paste0("assets/DESeq2results-", fileName, "_fc.tsv")
  write_delim(outTable, file=newName, delim = "\t")
  upName <- paste0(fileName, ".up")
  upName <-outTable %>% filter(log2FoldChange >= 1) %>%
  filter(padj <= 0.05) %>%
  arrange(desc(log2FoldChange)) %>%
  dplyr::select(Gene)
  downName <- paste0(fileName, ".donw")
  downName <-outTable %>% filter(log2FoldChange <= 1) %>%
  filter(padj <= 0.05) %>%
  arrange(desc(log2FoldChange)) %>%
  dplyr::select(Gene)
}

The results are saved as tsv files

processDE(res.nCTvsBAP, "nCTvsBAP")
processDE(res.nCTvsSTB, "nCTvsSTB")
processDE(res.nTEvsSTB, "nTEvsSTB")
processDE(res.nTEvsBAP, "nTEvsBAP")
processDE(res.nTEvsL40, "nTEvsL40")
processDE(res.nCTvsL40, "nCTvsL40")
processDE(res.BAPvsL40, "BAPvsL40")
processDE(res.STBvsL40, "STBvsL40")
processDE(res.BAPvsSTB, "BAPvsSTB")

Preparing GeneLists for PCE

We will create a function to save the DE gene lists (up and down) and run PCE on them we will need the annotations to convert the gene-id-version to just gene-id (and add other metadata later)

ensembl=useMart("ENSEMBL_MART_ENSEMBL")
listDatasets(ensembl) %>%
  filter(str_detect(description, "Human"))
#>                 dataset              description    version
#> 1 hsapiens_gene_ensembl Human genes (GRCh38.p13) GRCh38.p13
ensembl = useDataset("hsapiens_gene_ensembl", mart=ensembl)
listFilters(ensembl) %>%
  filter(str_detect(name, "ensembl"))
#>                            name
#> 1               ensembl_gene_id
#> 2       ensembl_gene_id_version
#> 3         ensembl_transcript_id
#> 4 ensembl_transcript_id_version
#> 5            ensembl_peptide_id
#> 6    ensembl_peptide_id_version
#> 7               ensembl_exon_id
#>                                                      description
#> 1                       Gene stable ID(s) [e.g. ENSG00000000003]
#> 2       Gene stable ID(s) with version [e.g. ENSG00000000003.15]
#> 3                 Transcript stable ID(s) [e.g. ENST00000000233]
#> 4 Transcript stable ID(s) with version [e.g. ENST00000000233.10]
#> 5                    Protein stable ID(s) [e.g. ENSP00000000233]
#> 6     Protein stable ID(s) with version [e.g. ENSP00000000233.5]
#> 7                              Exon ID(s) [e.g. ENSE00000000003]
filterType <- "ensembl_gene_id_version"
filterValues <- rownames(cts)
listAttributes(ensembl) %>%
  head(20)
#>                             name                                description
#> 1                ensembl_gene_id                             Gene stable ID
#> 2        ensembl_gene_id_version                     Gene stable ID version
#> 3          ensembl_transcript_id                       Transcript stable ID
#> 4  ensembl_transcript_id_version               Transcript stable ID version
#> 5             ensembl_peptide_id                          Protein stable ID
#> 6     ensembl_peptide_id_version                  Protein stable ID version
#> 7                ensembl_exon_id                             Exon stable ID
#> 8                    description                           Gene description
#> 9                chromosome_name                   Chromosome/scaffold name
#> 10                start_position                            Gene start (bp)
#> 11                  end_position                              Gene end (bp)
#> 12                        strand                                     Strand
#> 13                          band                             Karyotype band
#> 14              transcript_start                      Transcript start (bp)
#> 15                transcript_end                        Transcript end (bp)
#> 16      transcription_start_site             Transcription start site (TSS)
#> 17             transcript_length Transcript length (including UTRs and CDS)
#> 18                transcript_tsl             Transcript support level (TSL)
#> 19      transcript_gencode_basic                   GENCODE basic annotation
#> 20             transcript_appris                          APPRIS annotation
#>            page
#> 1  feature_page
#> 2  feature_page
#> 3  feature_page
#> 4  feature_page
#> 5  feature_page
#> 6  feature_page
#> 7  feature_page
#> 8  feature_page
#> 9  feature_page
#> 10 feature_page
#> 11 feature_page
#> 12 feature_page
#> 13 feature_page
#> 14 feature_page
#> 15 feature_page
#> 16 feature_page
#> 17 feature_page
#> 18 feature_page
#> 19 feature_page
#> 20 feature_page
attributeNames <- c('ensembl_gene_id_version', 
                    'ensembl_gene_id', 
                    'external_gene_name')
annot <- getBM(attributes=attributeNames,
               filters = filterType,
               values = filterValues,
               mart = ensembl)
# remove duplicates, if any
isDup <- duplicated(annot$ensembl_gene_id)
dup <- annot$ensembl_gene_id[isDup]
annot <- annot[!annot$ensembl_gene_id%in%dup,]

function to create up regulated gene list

pceUP <- function(restable) {
  restable <- restable[order(restable$padj), ]
  outTable <- merge(as.data.frame(restable), 
                    as.data.frame(counts(dds, normalized=TRUE)),
                    by="row.names", sort=FALSE)
  names(outTable)[1] <- "Gene"
  upName <-outTable %>% filter(log2FoldChange >= 1) %>%
  filter(padj <= 0.05) %>%
  arrange(desc(log2FoldChange)) %>%
  dplyr::select(Gene)
  upNew <- annot[annot$ensembl_gene_id_version %in% upName$Gene,]
  upList <- upNew$ensembl_gene_id
  return(upList)
}

function to create down regulated gene list

pceDOWN <- function(restable) {
  restable <- restable[order(restable$padj), ]
  outTable <- merge(as.data.frame(restable),
                    as.data.frame(counts(dds, normalized=TRUE)), 
                    by="row.names", sort=FALSE)
  names(outTable)[1] <- "Gene"
  downName <-outTable %>% filter(log2FoldChange <= -1) %>%
  filter(padj <= 0.05) %>%
  arrange(desc(log2FoldChange)) %>%
  dplyr::select(Gene)
  downNew <- annot[annot$ensembl_gene_id_version %in% downName$Gene,]
  downList <- downNew$ensembl_gene_id
  return(downList)
}

Run the functions on the DESeq2 results object

nCTvsBAP.up <- pceUP(res.nCTvsBAP)
nCTvsSTB.up <- pceUP(res.nCTvsSTB)
nTEvsSTB.up <- pceUP(res.nTEvsSTB)
nTEvsBAP.up <- pceUP(res.nTEvsBAP)
nTEvsL40.up <- pceUP(res.nTEvsL40)
nCTvsL40.up <- pceUP(res.nCTvsL40)
BAPvsL40.up <- pceUP(res.BAPvsL40)
STBvsL40.up <- pceUP(res.STBvsL40)
BAPvsSTB.up <- pceUP(res.BAPvsSTB)
nCTvsBAP.down <- pceDOWN(res.nCTvsBAP)
nCTvsSTB.down <- pceDOWN(res.nCTvsSTB)
nTEvsSTB.down <- pceDOWN(res.nTEvsSTB)
nTEvsBAP.down <- pceDOWN(res.nTEvsBAP)
nTEvsL40.down <- pceDOWN(res.nTEvsL40)
nCTvsL40.down <- pceDOWN(res.nCTvsL40)
BAPvsL40.down <- pceDOWN(res.BAPvsL40)
STBvsL40.down <- pceDOWN(res.STBvsL40)
BAPvsSTB.down <- pceDOWN(res.BAPvsSTB)

Run PCE

The above gene lists are used for running PCE. The function used for running PCE is below

# load the PCE data
l <- load(file = "~/TutejaLab/PlacentaEnrich/combine-test-expression1.Rdata")
humanGeneMapping <- dataset$GRCH38$humanGeneMapping
d <- dataset$PlacentaDeciduaBloodData
data <- d$expressionData
cellDetails <- d$cellDetails
# create a run PCE function
runpce <- function(inputgenelist, title, shade) {
  inputGenes<-toupper(inputgenelist)
  expressionData<-data[intersect(row.names(data),humanGeneMapping$Gene),]
  se<-SummarizedExperiment(assays = SimpleList(as.matrix(expressionData)),
                           rowData = row.names(expressionData),
                           colData = colnames(expressionData))
  cellSpecificGenesExp<-teGeneRetrieval(se,expressedGeneThreshold = 1)
  print(length(inputGenes))
  gs<-GeneSet(geneIds=toupper(inputGenes))
  output2<-teEnrichmentCustom(gs,cellSpecificGenesExp)
  enrichmentOutput<-setNames(data.frame(assay(output2[[1]]),
                                        row.names = rowData(output2[[1]])[,1]),
                             colData(output2[[1]])[,1])
  row.names(cellDetails)<-cellDetails$RName
  enrichmentOutput$Tissue<- cellDetails[row.names(enrichmentOutput),"CellName"]
  ggplot(data = enrichmentOutput, mapping = aes(x = reorder (Tissue, -Log10PValue), Log10PValue)) +
    geom_bar(stat = "identity", color = shade, fill = shade) + theme_classic(base_size = 10) + 
    theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1, size = 10),
          plot.title = element_text(size = 12),
          plot.margin = unit(c(1,1,1,2), "cm")) +
          labs(x = "",
               y = "-log10 p-value (adj.)",
               title = title) +
    scale_y_continuous(expand = expansion(mult = c(0, .1)))
}

The PCE is run on each of the gene list as follows (up and down pair are displayed together)

p <- runpce(nCTvsBAP.up , "overexpressed in H9_nCT_P10_Io", "#F16746")
#> [1] 1381
#> No background list provided. Using all the
#>                 genes as background.
q <- runpce(nCTvsBAP.down , "overexpressed in pBAP_D3_Io", "#0773B2")
#> [1] 2821
#> No background list provided. Using all the
#>                 genes as background.
panel_plot <- plot_grid(p, q, labels=c("X", "Y"), ncol=1, nrow = 2)
panel_plot

FigX: H9_nCT_P10_Io vs. BAP.d3

p <- runpce(nCTvsSTB.up , "overexpressed in H9_nCT_P10_Io)", "#F16746")
#> [1] 563
#> No background list provided. Using all the
#>                 genes as background.
q <- runpce(nCTvsSTB.down , "overexpressed in H1_BAP_D8_>70_Yabe", "#5C823A")
#> [1] 1766
#> No background list provided. Using all the
#>                 genes as background.
panel_plot <- plot_grid(p, q, labels=c("X", "Y"), ncol=1, nrow = 2)
panel_plot

FigX: H9_nCT_P10_Io vs. H1_BAP_D8_>70_Yabe

p <- runpce(nTEvsSTB.up , "overexpressed in H9_nTE_D3_Io", "#C79D2E")
#> [1] 1069
#> No background list provided. Using all the
#>                 genes as background.
q <- runpce(nTEvsSTB.down , "overexpressed in H1_BAP_D8_>70_Yabe", "#5C823A")
#> [1] 2070
#> No background list provided. Using all the
#>                 genes as background.
panel_plot <- plot_grid(p, q, labels=c("X", "Y"), ncol=1, nrow = 2)
panel_plot

FigX: H9_nTE_D3_Io vs. H1_BAP_D8_>70_Yabe

p <- runpce(nTEvsBAP.up , "overexpressed in H9_nTE_D3_Io", "#C79D2E")
#> [1] 886
#> No background list provided. Using all the
#>                 genes as background.
q <- runpce(nTEvsBAP.down , "overexpressed in pBAP_D3_Io", "#0773B2")
#> [1] 2160
#> No background list provided. Using all the
#>                 genes as background.
panel_plot <- plot_grid(p, q, labels=c("X", "Y"), ncol=1, nrow = 2)
panel_plot

FigX: H9_nTE_D3_Io vs. pBAP_D3_Io

p <- runpce(nTEvsL40.up , "overexpressed in H9_nTE_D3_Io", "#C79D2E")
#> [1] 1183
#> No background list provided. Using all the
#>                 genes as background.
q <- runpce(nTEvsL40.down , "overexpressed in H1_BAP_D8_<40_Yabe", "#AFBF38")
#> [1] 2260
#> No background list provided. Using all the
#>                 genes as background.
panel_plot <- plot_grid(p, q, labels=c("X", "Y"), ncol=1, nrow = 2)
panel_plot

FigX: H9_nTE_D3_Io vs. H1_BAP_D8_<40_Yabe

p <- runpce(nCTvsL40.up , "overexpressed in H9_nCT_P10_Io", "#F16746")
#> [1] 650
#> No background list provided. Using all the
#>                 genes as background.
q <- runpce(nCTvsL40.down , "overexpressed in H1_BAP_D8_<40_Yabe", "#AFBF38")
#> [1] 2095
#> No background list provided. Using all the
#>                 genes as background.
panel_plot <- plot_grid(p, q, labels=c("X", "Y"), ncol=1, nrow = 2)
panel_plot

FigX: H9_nCT_P10_Io vs. H1_BAP_D8_<40_Yabe

p <- runpce(BAPvsL40.up , "up regulated in pBAP_D3_Io", "#0773B2")
#> [1] 1544
#> No background list provided. Using all the
#>                 genes as background.
q <- runpce(BAPvsL40.down , "overexpressed in H1_BAP_D8_<40_Yabe", "#AFBF38")
#> [1] 1351
#> No background list provided. Using all the
#>                 genes as background.
panel_plot <- plot_grid(p, q, labels=c("X", "Y"), ncol=1, nrow = 2)
panel_plot

FigX: pBAP_D3_Io vs. H1_BAP_D8_<40_Yabe

p <- runpce(STBvsL40.up , "overexpressed in H1_BAP_D8_>70_Yabe", "#AFBF38")
#> [1] 1266
#> No background list provided. Using all the
#>                 genes as background.
q <- runpce(STBvsL40.down , "overexpressed in H1_BAP_D8_<40_Yabe", "#0773B2")
#> [1] 1157
#> No background list provided. Using all the
#>                 genes as background.
panel_plot <- plot_grid(p, q, labels=c("X", "Y"), ncol=1, nrow = 2)
panel_plot

FigX: H1_BAP_D8_>70_Yabe vs. H1_BAP_D8_<40_Yabe

p <- runpce(BAPvsSTB.up , "overexpressed in pBAP_D3_Io", "#0773B2")
#> [1] 1732
#> No background list provided. Using all the
#>                 genes as background.
q <- runpce(BAPvsSTB.down , "overexpressed in H1_BAP_D8_>70_Yabe", "#5C823A")
#> [1] 1424
#> No background list provided. Using all the
#>                 genes as background.
panel_plot <- plot_grid(p, q, labels=c("X", "Y"), ncol=1, nrow = 2)
panel_plot

FigX: pBAP_D3_Io vs. H1_BAP_D8_>70_Yabe

Volcano plots

We will setup a funciton for drawing volcano plots and then run them on the DESeq2 results. The funciton is as shown below

mart <- read.csv("~/TutejaLab/amnion_for_ms_20210715/de-analyses/mart-genes.tsv",
                 sep="\t", stringsAsFactors = TRUE,
                 header = TRUE)
myVolPlot <- function(restable, first, second, shade1, shade2) {
  restable <- restable[order(restable$padj), ]
  outTable <- tibble::rownames_to_column(as.data.frame(restable), "Gene")
  outTable <- merge(outTable, mart, by.x=c("Gene"), 
                    by.y = c("ensembl_gene_id_version"),
                    sort=FALSE)
  outTable <- outTable %>% 
    mutate_all(na_if,"") %>% 
    mutate_all(na_if," ") %>% 
    mutate(gene_symbol = coalesce(gene_symbol,Gene))
  outTable$diffexpressed <- "other.genes"
  outTable$diffexpressed[outTable$log2FoldChange >= 1 &
                           outTable$padj <= 0.05] <- paste0("Higher expression in ", first)
  outTable$diffexpressed[outTable$log2FoldChange <= -1 &
                           outTable$padj <= 0.05] <- paste0("Higher expression in ", second)
  outTable$delabel <- ""
  outTable$delabel[outTable$log2FoldChange >= 1 
             & outTable$padj <= 0.05 
             & !is.na(outTable$padj)] <- outTable$gene_symbol[outTable$log2FoldChange >= 1 
                                                      & outTable$padj <= 0.05 
                                                      & !is.na(outTable$padj)]
  outTable$delabel[outTable$log2FoldChange <= -1 
             & outTable$padj <= 0.05 
             & !is.na(outTable$padj)] <- outTable$gene_symbol[outTable$log2FoldChange <= -1 
                                                      & outTable$padj <= 0.05 
                                                      & !is.na(outTable$padj)]
  ggplot(outTable, aes(x=log2FoldChange, y=-log10(padj), col=diffexpressed, label=delabel)) +
  geom_point(alpha = 0.5) +
  theme_classic() +
  scale_color_manual(name = "Expression", values=c(shade1, shade2, "#4d4d4d")) +
#  ggtitle(paste0(first, " vs. ", second)) +
  xlab("log2 fold change") +
  ylab("-log10 pvalue (adjusted)") +
  theme(legend.text.align = 0,
    legend.position = c(.95, .95),
    legend.justification = c("right", "top"),
    legend.box.just = "right",
    legend.margin = margin(5, 5, 5, 5))
}

Running Volcano plots for each comparison is shown below

ggplotly(myVolPlot(res.nCTvsBAP, "H9_nCT_P10_Io", "H9_pBAP_D3_Io", "#F16746", "#0571b0" ))

FigX: H9_nCT_P10_Io vs. H9_pBAP_D3_Io

ggplotly(myVolPlot(res.nCTvsSTB, "H9_nCT_P10_Io", "H1_BAP_D8_>70_Yabe", "#5C823A", "#F16746" ))

FigX: H9_nCT_P10_Io vs. H1_BAP_D8_>70_Yabe

ggplotly(myVolPlot(res.nTEvsSTB, "H9_nTE_D3_Io", "H1_BAP_D8_>70_Yabe", "#5C823A", "#C79D2E" ))

FigX: H9_nTE_D3_Io vs. H1_BAP_D8_>70_Yabe

ggplotly(myVolPlot(res.nTEvsBAP, "H9_nTE_D3_Io", "H9_pBAP_D3_Io", "#C79D2E", "#0571b0" ))

FigX: H9_nTE_D3_Io vs. H9_pBAP_D3_Io

ggplotly(myVolPlot(res.nTEvsL40, "H9_nTE_D3_Io", "H1_BAP_D8_<40_Yabe", "#AEBD38", "#C79D2E"))

FigX: H9_nTE_D3_Io vs. H1_BAP_D8_<40_Yabe

ggplotly(myVolPlot(res.nCTvsL40, "H9_nCT_P10_Io", "H1_BAP_D8_<40_Yabe", "#AEBD38", "#F16746" ))

FigX: H9_nCT_P10_Io vs. H1_BAP_D8_<40_Yabe

ggplotly(myVolPlot(res.BAPvsL40, "H9_pBAP_D3_Io", "H1_BAP_D8_<40_Yabe", "#AEBD38", "#0571b0" ))

FigX: H9_pBAP_D3_Io vs. H1_BAP_D8_<40_Yabe

ggplotly(myVolPlot(res.STBvsL40, "H1_BAP_D8_>70_Yabe", "H1_BAP_D8_<40_Yabe", "#AEBD38", "#5C823A" ))

FigX: H1_BAP_D8_>70_Yabe vs. H1_BAP_D8_<40_Yabe

ggplotly(myVolPlot(res.BAPvsSTB, "H9_pBAP_D3_Io", "H1_BAP_D8_>70_Yabe", "#5C823A", "#0571b0" ))

FigX: H9_pBAP_D3_Io vs. H1_BAP_D8_>70_Yabe

Section 2: DESeq2 analysis

Arun Seetharam, Ha Vu

2021-08-13